{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "if IN_COLAB:\n",
    "    !git clone https://github.com/cs357/demos-cs357.git\n",
    "    !mv demos-cs357/figures/ .\n",
    "    !mv demos-cs357/additional_files/ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving linear system of equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Useful libraries for this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import scipy.linalg as sla\n",
    "import scipy.sparse as sparse\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from PIL import Image\n",
    "\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=2)\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "\n",
    "from time import time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Transforming images using linear operators:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import an image of a Social Security Number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "img = Image.open('figures/ssn.png')\n",
    "xmat = (255 - np.asarray(img).max(axis=2))/255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(xmat.shape)\n",
    "print(xmat.min(),xmat.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.imshow(xmat);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \"Vectorize\" the image, creating the 1d array `x`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = xmat.flatten()\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct a \"blur\" matrix\n",
    "More about this blur matrix on a later MP..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "imat, jmat = np.meshgrid(np.arange(xmat.shape[0]), np.arange(xmat.shape[1]), indexing='ij')\n",
    "\n",
    "ivec = np.atleast_2d(imat.flatten())\n",
    "jvec = np.atleast_2d(jmat.flatten())\n",
    "\n",
    "A = np.fmax(0, 1 - np.sqrt((ivec.T - ivec)**2 + (jvec.T - jvec)**2)/5)\n",
    "A /= A.sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute b = A x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = A @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b2D=b.reshape(xmat.shape)\n",
    "plt.imshow(b2D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assume we have the blurred image (b), solve for the unblurred one (x)\n",
    "\n",
    "In essence, we want $x = A^{-1} b$. Instead of evaluating the inverse, we can evaluate:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### a)  `np.linalg.solve`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_solve1 = la.solve(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.imshow(x_solve1.reshape(xmat.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila! We have our recovered original image :-)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What is happening inside \"solve\"?\n",
    "\n",
    "A factorization and substitutions. For example, let's think of LU factorization:\n",
    "\n",
    "If $Ax = P L U x = b$, then there are two steps:\n",
    "1. $y \\leftarrow \\text{solve}\\,\\, L y = P^Tb$\n",
    "2. $x \\leftarrow \\text{solve}\\,\\, U x = y$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### b) LU factorization and triangular substitutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.linalg as sla\n",
    "P, L, U = sla.lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# L y = P^T b\n",
    "y = sla.solve_triangular(L, P.T@b, lower=True)\n",
    "\n",
    "# U x = y\n",
    "x_solve = sla.solve_triangular(U, y)\n",
    "\n",
    "plt.imshow(x_solve.reshape(xmat.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Why not just `np.linalg.solve`?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you have many social security numbers to un-blur. You should factorize your blur matrix only once and then perform several triangular solves.\n",
    "\n",
    "Let's time:\n",
    "\n",
    "1. solve including factorization\n",
    "2. factorization\n",
    "3. solve, given a factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit P, L, U = sla.lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.solve_triangular(U, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's take a look at the matrix format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(131)\n",
    "plt.spy(A); plt.xticks([]); plt.yticks([]);\n",
    "\n",
    "plt.subplot(132)\n",
    "plt.spy(L); plt.xticks([]); plt.yticks([]);\n",
    "\n",
    "plt.subplot(133)\n",
    "plt.spy(U); plt.xticks([]); plt.yticks([]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that many entries of the matrices are zeros! We don't need to store all of the zeros. In fact, we can even use solvers that efficiently take into account the format of the matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    " A_csr = sparse.csr_matrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sparse.linalg.spsolve(A_csr,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_solve_3 = sparse.linalg.spsolve(A_csr,b)\n",
    "plt.imshow(x_solve_3.reshape(xmat.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is just a teaser :-) We will talk more later about sparse systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now add some noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b_noisy = b + 1e-5 * np.random.rand(b.size)\n",
    "\n",
    "y = sla.solve_triangular(L, np.dot(P.T, b_noisy), lower=True)\n",
    "x_solve = sla.solve_triangular(U, y)\n",
    "\n",
    "plt.imshow(x_solve.reshape(xmat.shape))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "What happened to your solution. We will soon be talking about conditioning of matrices. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  2) Relative cost of matrix operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_values = (10**np.linspace(1, 4, 15)).astype(np.int32)\n",
    "n_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "times_matmul = []\n",
    "times_lu = []\n",
    "\n",
    "for n in n_values:\n",
    "    print(n)\n",
    "    A = np.random.randn(n, n)\n",
    "    start_time = time()\n",
    "    A.dot(A)\n",
    "    times_matmul.append(time() - start_time)\n",
    "    start_time = time()\n",
    "    sla.lu(A)\n",
    "    times_lu.append(time() - start_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(n_values, times_matmul, label='matmul')\n",
    "plt.plot(n_values, times_lu, label='lu')\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"Matrix size $n$\")\n",
    "plt.ylabel(\"Wall time [s]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The faster algorithms make the slower ones look bad. But... it's all relative.\n",
    "* Can we get a better plot?\n",
    "* Can we see the asymptotic cost ($O(n^3)$) of these algorithms from the plot?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.loglog(n_values, times_matmul, label='matmul')\n",
    "plt.loglog(n_values, times_lu, label='lu')\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"Matrix size $n$\")\n",
    "plt.ylabel(\"Wall time [s]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Example:  Topology design optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](figures\\myfigure.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When performing optimization of structural problem, for example to obtain the bridge design above, you may want to use a numerical method called Finite Element Method. The optimization will consist of a series of `solve` of the form:\n",
    "$$ {\\bf K} {\\bf u} = {\\bf F} $$\n",
    "\n",
    "Here will load the matrix $ {\\bf K}$ from a file. The matrix is given in Compressed Sparse Column (CSC) format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "K = sparse.load_npz('additional_files/yourmatrix.npz')\n",
    "K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can `spy` the distribution of the non-zero entries of the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.spy(K)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix ${\\bf K}$ has a banded format, and it is also symmetric and positive definite. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Kdense = K.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Kdense.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.max(Kdense-Kdense.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sla.norm(Kdense-Kdense.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the linear system of equations using different methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F = np.zeros(K.shape[0])\n",
    "F[1]=-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### a) la.solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u1 = sla.solve(Kdense,F)\n",
    "u1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.solve(Kdense,F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### b) LU factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "P,L,U = sla.lu(Kdense)\n",
    "y = sla.solve_triangular(L, np.dot(P.T, F), lower=True)\n",
    "u2 = sla.solve_triangular(U, y)\n",
    "u2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.lu(Kdense)\n",
    "%timeit sla.solve_triangular(L, np.dot(P.T, F), lower=True)\n",
    "%timeit sla.solve_triangular(U, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### c) Cholesky factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Kcho = sla.cholesky(Kdense)\n",
    "u3 = sla.cho_solve((Kcho,False),F)\n",
    "u3.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.cholesky(Kdense)\n",
    "%timeit sla.cho_solve((Kcho,False),F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### d) Sparse solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import spsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u4 = spsolve(K,F)\n",
    "u4.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit spsolve(K,F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Steady-State Advection Diffusion Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('./additional_files')\n",
    "from CA_4_support import *\n",
    "from matplotlib import animation,rc\n",
    "from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Consider the following ordinary differential equation (ODE) for $x$ in the interval $[-1,1]$:\n",
    "\n",
    "$$-u''(x) + u'(x) = f(x)$$\n",
    "\n",
    "With boundary conditions\n",
    "\n",
    "$$u(-1) = u(1) = 0$$\n",
    "           \n",
    "Here, $f(x)$ is a known function and $u(x)$ is the unknown that we are trying to solve for."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral Element Approximation\n",
    "The spectral element method (SEM) is a way to solve such an ODE numerically.  In the simplest case, SEM approximates $u(x)$ by a polynomial of degree $N$.\n",
    "\n",
    "A polynomial can be uniquely described by $P = N+1$ points on the interval $[-1,1]$.  So the SEM introduces a set of $P$ discrete points\n",
    "\n",
    "$-1 = x_0 < x_1 < x_2 < \\dots < x_{N-1} < x_N = 1$\n",
    "\n",
    "and constructs a linear system\n",
    "\n",
    "$\\mathbf{A}\\mathbf{u} = \\mathbf{b}$\n",
    "\n",
    "that is solved to obtain the value of $u$ at each point $x_i$. Here, the vector $\\mathbf{u} = [u_1,u_2,\\dots,u_{N-1}]^T$, where $u_i = u(x_i)$ for $i = 1,2,\\dots,N-1$.\n",
    "\n",
    "Note that we exclude the points $x_0=-1$ and $x_N = 1$ because we already know the value of $u_0$ and $u_N$ from the boundary conditions.\n",
    "\n",
    "$\\mathbf{b}$ is a vector that depends on the function $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll provide a function that takes a polynomial order $N$ as input and creates both the partition points $x_i$ defined in the interval $([-1,1]$, and the system matrix $\\mathbf{A}$.  We'll also provide a function that takes the  points $x_i$ and a known function $f(x)$ and returns the vector $\\mathbf{b}$. The function signatures are defined below:\n",
    "\n",
    "```python\n",
    "A,x =  SEM_system_1(N)\n",
    "\n",
    "def f(z):\n",
    "    # define the function f as a function of z\n",
    "    return f\n",
    "  \n",
    "b = SEM_rhs_1(f,x)\n",
    "```\n",
    "Generate the arrays $\\mathbf{A}$, $\\mathbf{b}$, and $\\mathbf{x}$ using the provided functions. Use for example `N = 30`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 30  # polynomial order\n",
    "A,x = SEM_system_1(N)\n",
    "\n",
    "def f(z):\n",
    "    return 2*(z-1)\n",
    "\n",
    "b = SEM_rhs_1(f,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check out the shapes of all 3 objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array of points $\\mathbf{x}$ is of size $N+1$ but the matrix and right hand side lead to a linear system of size $(N-1) \\times (N-1)$.  This is because we don't need to solve for the value of $u$ at the two endpoints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use $\\mathbf{A}$, $\\mathbf{b}$, and $\\mathbf{x}$ to solve for $\\mathbf{u}$.  Plot your numerical approximation and the exact solution $u_e(x) = x^2 - 1$.  Don't forget to add zeros to the start and end of the vector to account for the boundary conditions.  Do the two solutions match up?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "u = ...\n",
    "u_plot = ...\n",
    "u_exact = x**2 - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "plt.plot(x,u_plot,'o')\n",
    "plt.plot(x,u_exact,'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run the cells above changing the number of points. Use larger and smaller values. What do you notice?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution will change for different right hand sides.  Experiment with different functions $f(x)$ and see what the solution looks like. For example, try the function:\n",
    "\n",
    "$$ f(x) = \\sin(1.5 \\pi x) $$\n",
    "\n",
    "Note that since your ODE has not changed, your system matrix  $\\mathbf{A}$ is still the same! You only need to redefine the vector $\\mathbf{b}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(z):\n",
    "    return ...\n",
    "\n",
    "b = SEM_rhs_1(f,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Obtain u and plot \n",
    "u =...\n",
    "u_plot = ...\n",
    "\n",
    "plt.plot(x,u_plot,'-o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5) Time-Dependent Advection Diffusion\n",
    "\n",
    "Now we'll look at the time-dependent advection diffusion equation, which is a partial differential equation (PDE)\n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t} + \\frac{\\partial u}{\\partial x} - D\\frac{\\partial^2 u}{\\partial x^2} = 0$$\n",
    "\n",
    "with boundary conditions:\n",
    "\n",
    "$$u(-1,t) = 0 = \\frac{\\partial u}{\\partial x}(1,t)$$\n",
    "\n",
    "Where $D = 0.1$ is the diffusion coefficient.  Notice how the boundary condition at the right hand side is a little different than the first equation.\n",
    "\n",
    "In order to solve this equation we again use the SEM.  However, the equation is time dependent, so we also partition the time dimension:\n",
    "\n",
    "$$0 = t_0 < t_1 < t_2 < \\dots$$.\n",
    "\n",
    "We assume the $t_i$'s are evenly spaced with distance $\\Delta t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the solution $u$ depends on both space and time: $u(x,t)$.  For each time $t_n$, we can define a vector of values of $u$ like above:\n",
    "\n",
    "$\\mathbf{u}^{(n)} = [u_1^{(n)},\\dots, u_N^{(n)}]^T$.  Here $u_i^{(n)} = u(x_i,t_n)$.  We assume we know $u(x,0)$ so that $\\mathbf{u}^{(0)}$ is known.  Then we get a *timestepping scheme*:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{A}\\mathbf{u}^{(1)} &= \\mathbf{b}^{(0)}\\\\\n",
    "\\mathbf{A}\\mathbf{u}^{(2)} &= \\mathbf{b}^{(1)}\\\\\n",
    "\\mathbf{A}\\mathbf{u}^{(3)} &= \\mathbf{b}^{(2)}\\\\\n",
    "&\\ \\vdots\\\\\n",
    "\\mathbf{A}\\mathbf{u}^{(n+1)} &= \\mathbf{b}^{(n)}\n",
    "\\end{align}\n",
    "\n",
    "Here $\\mathbf{b}^{(n)}$ is a vector that depends on $\\mathbf{u}^{(n)}$, and can be obtained using the helper function:\n",
    "```\n",
    "b = SEM_rhs_2(un)\n",
    "```\n",
    "We'll select a value of $\\Delta t = 0.01$ and take 200 timesteps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dt = 0.01        # time increment\n",
    "timesteps = 200  # number of time steps\n",
    "N = 50           # polynomial order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix $\\mathbf{A}$ depends not just on the polynomial order $N$ but also the value of $\\Delta t$. We provide the helper function:\n",
    "\n",
    "```python\n",
    "A,x = SEM_system_2(N,dt)\n",
    "```\n",
    "Generate the arrays $\\mathbf{A}$, and $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A,x = SEM_system_2(N,dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also provide a function that defines the initial condition $u(x,0)$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def u_initial(z):\n",
    "    a = -10*(z+1)*z**6*(z-1)*(z< 0)\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will store all the vectors $\\mathbf{u}^{(0)}$, $\\mathbf{u}^{(1)}$, $\\mathbf{u}^{(2)}$,$\\dots$ as columns of a matrix. Initialize with zeros the 2d numpy array (matrix) $\\mathbf{u}$ with the appropriate shape. Update the first column with the initial condition using the function `u_initial`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u = np.zeros((x.shape[0],timesteps+1))\n",
    "# store the initial condition in the zeroth column\n",
    "u[:,0] = u_initial(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# plot initial condition\n",
    "plt.plot(x,u[:,0],'o-')\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$u$')\n",
    "plt.title('Initial Condition')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the array $\\mathbf{b}$ corresponding to the initial condition `u[:,0]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First time step\n",
    "u_old = u[:,0]\n",
    "b = SEM_rhs_2(u_old)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take the first time step and store the solution in `u[:,1]`. Make sure you are handling the size of your arrays correctly and including the boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u_new = np.linalg.solve(A,b)\n",
    "u[:,1] = np.append(0,u_new)   # add back boundary condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can take the second time step, and store your solution in `u[:,2]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# second time step\n",
    "u_old = u[:,1]\n",
    "b = SEM_rhs_2(u_old)\n",
    "u_new = np.linalg.solve(A,b)\n",
    "u[:,2] = np.append(0,u_new) # add back boundary condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot what you have so far:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(x,u[:,0],label = '$\\mathbf{u}^{(0)}$')\n",
    "plt.plot(x,u[:,1],label ='$\\mathbf{u}^{(1)}$' )\n",
    "plt.plot(x,u[:,2],label = '$\\mathbf{u}^{(2)}$')\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$u$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only 198 time steps to go...\n",
    "We'll obviously want to do this in a loop.  But calling `np.linalg.solve` every time step is too expensive.  We'll want to use the LU decomposition of $A$ to make computing the solution less expensive.  Here $\\mathbf{A} = \\mathbf{PLU}$, where $\\mathbf{P}^{-1} = \\mathbf{P}^T$ and $\\mathbf{L}$ and $\\mathbf{U}$ are lower and upper diagonal matrices\n",
    "\n",
    "The equation\n",
    "\n",
    "$\\mathbf{A}\\mathbf{u}^{(n+1)} = \\mathbf{b}^{(n)}$\n",
    "\n",
    "is transformed into \n",
    "\n",
    "$\\mathbf{PLUu}^{(n+1)} = \\mathbf{b}^{(n)}$\n",
    "\n",
    "so \n",
    "\n",
    "$\\mathbf{u}^{(n+1)} = \\mathbf{U}^{-1}\\mathbf{L}^{-1}\\mathbf{P}^T\\mathbf{b}^{(n)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the LU factorization of $\\mathbf{A}$ using `scipy.linalg.lu` and invert $\\mathbf{L}$ and $\\mathbf{U}$ using `scipy.linalg.solve_triangular`\n",
    "\n",
    "Check out the documentation for `scipy.linalg.solve_triangular` before you use them.  In particular, pay attention to the argument `lower`, and `unit_diagonal`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `plt.spy` to plot the non-zero pattern of each matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.spy(P)\n",
    "plt.figure()\n",
    "plt.spy(L)\n",
    "plt.figure()\n",
    "plt.spy(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the time history of the solution using a movie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots()\n",
    "ax.set_xlim((-1,1))\n",
    "ax.set_ylim((0,1.0))\n",
    "line, = ax.plot([],[],lw=2)\n",
    "\n",
    "def init():\n",
    "    line.set_data([],[])\n",
    "    return (line,)\n",
    "\n",
    "def animate(i):\n",
    "    y = u[:,i]\n",
    "    line.set_data(x,y)\n",
    "    return (line,)\n",
    "\n",
    "anim = animation.FuncAnimation(fig, animate, init_func=init, frames=200, interval=50, blit=True)\n",
    "rc('animation',html='html5')\n",
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6) Time-Dependent Diffusion Equation\n",
    "Now we'll look at the closely related time-dependent diffusion equation\n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t} - D\\frac{\\partial^2 u}{\\partial x^2} = 0$$\n",
    "\n",
    "with boundary conditions:\n",
    "\n",
    "$$u(-1,t) = 0 = u(1,t)$$\n",
    "\n",
    "We'll get another time-stepping scheme of the same form:\n",
    "\\begin{align}\n",
    "\\mathbf{A}\\mathbf{u}^{(1)} &= \\mathbf{b}^{(0)}\\\\\n",
    "\\mathbf{A}\\mathbf{u}^{(2)} &= \\mathbf{b}^{(1)}\\\\\n",
    "\\mathbf{A}\\mathbf{u}^{(3)} &= \\mathbf{b}^{(2)}\\\\\n",
    "&\\ \\vdots\\\\\n",
    "\\mathbf{A}\\mathbf{u}^{(n+1)} &= \\mathbf{b}^{(n)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will again define the number of time steps, the time increment, the polynomial degree and use a helper function to construct the matrix $\\mathbf{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dt = 0.01\n",
    "timesteps = 200\n",
    "N = 50\n",
    "A,x = SEM_system_3(N,dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diffusion equation leads to a *symmetric positive-definite* matrix.  This means that $\\mathbf{A} = \\mathbf{A}^T$ and the eigenvalues of $\\mathbf{A}$ are positive. (You'll learn about eigenvalues next week).\n",
    "\n",
    "- Check if the matrix is symmetric:\n",
    "- Check if the matrix has all positive eigenvalues (you can use np.linalg.eigvals(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# check for symmetry\n",
    "print(np.allclose(A,A.T))\n",
    "\n",
    "# check that smallest eigenvalue is positive\n",
    "eig_vals = np.linalg.eigvals(A)\n",
    "print(np.min(eig_vals) > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of an LU factorization, we'll use the Cholesky factorization: $\\mathbf{A} = \\mathbf{U}^T\\mathbf{U}$, where $\\mathbf{U}$ is an upper triangular matrix. This is a more efficient factorization method for symmetric positive-definite matrices) - about half the number of floating operations of LU factorization.\n",
    "\n",
    "Let's compare the computational time for both methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 200\n",
    "A,x = SEM_system_3(N,dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.cholesky(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%timeit sla.lu(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the time evolution we'll set `N = 50`, we have the following system matrix $\\mathbf{A}$ and initial condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 50\n",
    "A,x = SEM_system_3(N,dt)\n",
    "\n",
    "def u_initial(z):\n",
    "    a = -10*(z+1)*z**6*(z-1)*(z< 0)\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will store all the vectors $\\mathbf{u}^{(0)}$, $\\mathbf{u}^{(1)}$, $\\mathbf{u}^{(2)}$,$\\dots$ as columns of a matrix. Initialize with zeros the 2d numpy array (matrix) $\\mathbf{u}$ with the appropriate shape. Update the first column with the initial condition using the function `u_initial`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u = np.zeros((x.shape[0],timesteps+1))\n",
    "\n",
    "u[:,0] = u_initial(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `sla.cholesky` and `sla.solve_triangular` to evolve the solution forward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "U = sla.cholesky(A)\n",
    "\n",
    "for i in range(timesteps):\n",
    "    u_old = u[:,i]\n",
    "    b = SEM_rhs_3(u_old)\n",
    "    b1 = sla.solve_triangular(U.T,b,lower=True,unit_diagonal=False)\n",
    "    u_new = sla.solve_triangular(U,b1,lower=False,unit_diagonal=False)\n",
    "    u[:,i+1] = np.append(0,np.append(u_new,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another movie.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots()\n",
    "ax.set_xlim((-1,1))\n",
    "ax.set_ylim((0,1.0))\n",
    "line, = ax.plot([],[],lw=2)\n",
    "\n",
    "def init():\n",
    "    line.set_data([],[])\n",
    "    return (line,)\n",
    "\n",
    "def animate(i):\n",
    "    y = u[:,i]\n",
    "    line.set_data(x,y)\n",
    "    return (line,)\n",
    "\n",
    "anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                               frames=100, interval=50, blit=True)\n",
    "\n",
    "rc('animation',html='html5')\n",
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7) LU Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Factorize the matrix $A$ as $A = LU$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 4\n",
    "\n",
    "A = np.array([[2.0,8,4,1],[1,2,3,3],[1,2,6,2],[1,3,4,2]])\n",
    "\n",
    "## Algorithm 1\n",
    "## Factorization using the block-format, \n",
    "## creating new matrices L and U\n",
    "## and not modifying A\n",
    "print(\"LU factorization using Algorithm 1\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Algorithm 2\n",
    "## Factorization using the block-format \n",
    "## Matrices L and U are stored in the input matrix\n",
    "## that could be a copy of A or A itself\n",
    "print(\"LU factorization using Algorithm 2\")  \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a) Factorize and solve resulting triangular systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Factorize:\n",
    "\n",
    "$$A = L U$$\n",
    "\n",
    "to solve $Ax=b$\n",
    "\n",
    "2) Forward substitution:\n",
    "\n",
    "$y \\leftarrow \\text{solve}\\,\\, L y = b$\n",
    "    \n",
    "$$ y_i = b_i - \\sum_{j=1}^{i-1} L_{ij} y_j \\quad i=1,2,3,...,n$$\n",
    "\n",
    "Backward substitution:\n",
    "    \n",
    "$x \\leftarrow \\text{solve}\\,\\, U x = y$\n",
    "\n",
    "$$ x_i = y_i - \\sum_{j=i+1}^{n} U_{ij} x_j \\quad i=n, n-1, n-2, ..., 1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Factorize and solve \n",
    "M = A.copy()\n",
    "for i in range(n):\n",
    "    M[i+1:,i] = M[i+1:,i]/M[i,i]\n",
    "    M[i+1:,i+1:] -= np.outer(M[i+1:,i],M[i,i+1:]) \n",
    "\n",
    "print(\"Solve system of equations for \")\n",
    "b = np.array([2,2,1,4])\n",
    "print(\"b = \",b)\n",
    "x = np.zeros(n)\n",
    "y = np.zeros(n)\n",
    "\n",
    "# Foward-substitution - solve Ly = b for y\n",
    "for i in range(0, n):\n",
    "    tmp = b[i]\n",
    "    for j in range(0, i):\n",
    "        tmp -= M[i,j]*y[j]\n",
    "    y[i]=tmp\n",
    "#print(y)\n",
    "\n",
    "# Backward substitution - solve Ux = y for x\n",
    "for i in range(n-1, -1, -1):\n",
    "    tmp = y[i]\n",
    "    for j in range(i+1, n):\n",
    "        tmp -= M[i,j]*x[j]\n",
    "    x[i] = tmp/M[i,i]\n",
    "print(\"Solution is:\")\n",
    "print(\"x = \",x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b) Simple example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$$\n",
    "Ax = \\begin{bmatrix}c&1\\\\-1&1\\end{bmatrix}\\begin{bmatrix}x_1\\\\x_2\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}b_1\\\\b_2\\end{bmatrix}\n",
    "$$\n",
    "with an exact solution of\n",
    "$$\n",
    "x_{exact} = \\begin{bmatrix}1\\\\1\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def myLU(A): \n",
    "    M = A.copy()\n",
    "    for i in range(A.shape[0]):\n",
    "        M[i+1:,i] = M[i+1:,i]/M[i,i]\n",
    "        M[i+1:,i+1:] -= np.outer(M[i+1:,i],M[i,i+1:]) \n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Creating the arrays\n",
    "c = 1e-1\n",
    "A = np.array([[c, 1.], [-1, 1]])\n",
    "# xx is the exact solution\n",
    "xx = np.array([1,1])\n",
    "b = A.dot(xx)\n",
    "\n",
    "# Comput the LU\n",
    "M = myLU(A)\n",
    "L = np.tril(M,-1) + np.eye(2,2)\n",
    "U = np.triu(M)\n",
    "\n",
    "# Solve\n",
    "# x is the numerical (xhat)\n",
    "y = sla.solve_triangular(L, b, lower=True)\n",
    "x = sla.solve_triangular(U, y)\n",
    "\n",
    "print(\"Exact solution = \", xx)\n",
    "\n",
    "print(\"Computed solution = \",x)\n",
    "\n",
    "print(\"Error = \", la.norm(xx-x))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the matrix A singular?\n",
    "\n",
    "A) YES\n",
    "\n",
    "B) NO\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c) Pivoting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Creating the arrays\n",
    "c = 1e-1\n",
    "A = np.array([[c, 1.], [-1, 1]])\n",
    "# xx is the exact solution\n",
    "xx = np.array([1,1])\n",
    "b = A.dot(xx)\n",
    "\n",
    "\n",
    "# add code here\n",
    "\n",
    "print(\"Exact solution = \", xx)\n",
    "\n",
    "print(\"Computed solution = \",x)\n",
    "\n",
    "print(\"Error = \", la.norm(xx-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.linalg as sla\n",
    "import numpy as np\n",
    "\n",
    "n = 4\n",
    "\n",
    "A = np.array([[2.0,8,4,1],[1,2,3,3],[1,2,3,2],[1,3,4,2]])\n",
    "\n",
    "(P,L,U)=sla.lu(A)\n",
    "print(P)\n",
    "\n",
    "print(L)\n",
    "\n",
    "print(U)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
